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LONG  TERM  GOALS 

Development  of  accurate  and  fast  advanced  statistical  and  dynamical  nonlinear  models  of  ocean 
surface  waves,  based  on  first  physical  principles,  which  will  improve  and  accelerate  both  long  term 
ocean  surface  waves  forecasts  and  prediction  of  strongly  coherent  events,  such  as  freak  waves,  tsunami 
and  wave-breaking. 

OBJECTIVES 

Finding  of  physically  correct  wind  input  term  for  Hasselmann  equation,  understanding  of  the  balance 
of  source  terms  in  Hasselmann  equation,  investigation  of  the  problem  of  interaction  of  different  scales 
on  the  ocean  surface  (sea  and  swell),  development  of  new  water  surface  analytic  models  and  methods 
of  their  solution. 

APPROACH 

Advanced  analytical  techniques:  Hamiltonian  formalism,  self-similar  solutions,  analytic  solutions  of 
integral  equations;  numerical  methods  for  solution  of  integral  and  pseudo-differential  equations; 
comparison  of  analytic  and  numerical  results  with  experimental  data 

WORK  COMPLETED 

•  Confirmation  of  the  new  wind  input  term  through  the  limited  fetch  numerical  simulation 

•  Analytical  solution  and  its  numerical  confirmation  for  non-stationary  Hasselmann  equation 
without  non-linear  interaction  term 

•  Detection  of  swell  feedback  by  sea  background  through  theory  and  experimental  data 

•  Analytical  and  numerical  proof  of  non-integrability  of  2D  free-surface  hydrodynamics 

•  New  findings  in  mechanisms  of  the  energy  dissipation  during  wave  breaking  event 

RESULTS 

1.  Confirmation  of  the  new  wind  input  term  through  the  limited  fetch  numerical  simulation 
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The  Hasselmann  equation  for  energy  spectral  density  s(r,k,t) 


ds  doolr  ds 


dt  dk  dr 

in  the  limited  fetch  situation  can  be  reduced  to 


-  Snl  +  Sin  +  Sdiss 


^-co  s0^  =  S„+Swlm,+SJ„ 
ico  ox 


(2) 


where  x  is  the  coordinate  orthogonal  to  the  shore  and  9  is  the  angle  between  individual  wavenumber  k 
and  the  axis  x. 


Eq.(2)  is  somewhat  difficult  for  numerical  simulation,  as  it  contains  the  singularity  in  the  form  of 
(cos#)  on  the  right-hand  side.  We  overcame  this  problem  of  division  by  zero  through  zeroing  one 
half  of  the  Fourier  space  of  the  system  for  the  waves  propagating  toward  the  shore.  Since  it  is  well- 
known  that  the  energy  in  such  waves  is  small  with  respect  to  waves  propagating  in  the  offshore 
direction,  such  approximation  is  quite  reasonable  for  our  purposes. 

We  are  looking  for  wind  input  function  in  he  form  Swind  =  y(a>,  9)s(a>,  6)  where  y(co,  9)  =  acox+s  f  ( 9 ) . 


The  same  sort  of  self-similar  analysis  perfonned  for  the  time  domain  situation  (see  [1])  can  be  repeated 
for  the  limited  fetch  one.  The  result  of  inserting  the  self-similar  substitution 

s  =  xp+qF((oxq ) 

into  Eq.(2)  gives  the  values  of  the  indices 

\0q-\  1 

P  =  —^’  q  = 

2  s  2 

and 


q  = 


p  =  1,  s  = 


4 

3 


We  found  that  in  the  fetch-limited  case  the  wind  forcing  index  s  is  similar  to  the  time  domain  situation, 
and  the  wind  forcing  is  given  by  (see  [1]): 


S„w=0-2—  «>l 

P watt 


(  y/3 

CO 


V®o  j 


f(0),  f(0)  = 


Jcos  9  for  -  nil  <  6  <  nil 
[0  otherwise 


0)Q  =-M~,  -P™—  =  1.3  .  l(T3 

«10  P  water 


(3) 


To  check  the  theoretical  finding  (3),  the  numerical  simulation  of  the  Eq.(2)  has  been  perfonned  for 
Swind  input  function  (3),  dissipation  term  Sdiss  similar  to  [3],  where  white-capping  dissipation  tenn  was 

introduced  implicitly  through  f  5  (/  =  )  energy  spectral  tail  stretching  in  frequency  range  from 

2k 

fd  =1.1  to  /max  =  2.0 . 

Fig.  1  presents  characteristic  spectral  energy  distribution  for  this  limited  fetch  simulation. 
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Fig.l:  Energy  spectral  density  s{j\  6)  as  a  function  of  f  and  angle  0  in  polar  coordinates. 

Fig.2  presents  the  plot  of  the  function  fi  =  F(k )  •  k512  in  terms  of  the  parameter  (ujCp)1/3  /  g12  for  two 

different  runs,  corresponding  to  the  same  wind  speed  uw  =10.0  m/sec:  one  for  time  domain  simulation 

(see  [1]),  and  another  for  limited  fetch  simulation.  Both  simulations  show  good  correspondence  with 
Resio  et  al.  (2004)  [2]  regression  line. 
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Fig.2:  Two  simulations  for  wind  input  ulo=10.0  m/sec:  time-limited  domain  (triangles)  and  limited 

fetch  (crosses).  Dotted  line  -  correlation  of  the  equilibrium  range  coefficient  [)  with  (w2C/;  )13  /  g1  2 

based  on  data  from  six  disparate  sources  adopted  from  Resio  at  al,  2004  [2] .  Dashed  line  - 
theoretical  value  of  equilibrium  range  coefficient  () . 

2.  Analytical  solution  and  its  numerical  tests  for  Hasselmann  equation  without  nonlinear 
interaction 

The  Hasselmann  Eq.(l)  includes  essentially  different  physical  processes'  —  advection,  nonlinear 
interaction  and  external  forcing/damping.  It  is  important  to  understand  during  numerical 
implementation  of  Hasselmann  equation  simulation  that  individual  terms,  responsible  for  particular 
physical  process,  work  well.  In  this  connection  it's  necessary  to  test  each  of  them  individually  to  be 
sure  that  discretization  implementation  errors  do  not  overlap,  creating  wrongly  working  numerical 
code. 

As  far  as  concerns  nonlinear  interaction  term  Snl ,  we  can  be  sure  that  discrete  implementation  works 
quite  well.  Those  tests  include,  in  particular,  numerical  reproduction  of  analytical  Komogorov- 
Zakharov  and  self-similar  solutions  for  time-limited  fonnulation  of  Hasselmann  equation  without 
advection  term  [1],  [4]. 

The  question  of  testing  of  the  advection  part  of  Hasselmann  equation  (1)  remained  open  so  far.  We 
managed  to  find  specific  analytic  solution  and  show  numerically  that  it  is  reproduced  by  the 
corresponding  numerical  code. 

We  are  starting  with  linear  version  of  Hasselmann  Eq.(l)  stripped  off  the  nonlinear  interaction  term 
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ds  dco,  ds 
-  h - - - -ys 


dt  ok  dx 

For  one-dimensional  situation  and  deep  water  dispersion  relation^.  =  -Jgk  it  transforms  into 

ds  g  nds 

- 1 — cos  u —  =  ys 

dt  2  co  dx 


which  after  coordinate  transfonnation  X  — 


2  co 


g  cost? 


-x  becomes 


ds  ds 

- 1 - =  ys 

dt  dX 


\eytA<Xn 


Eq.(5)  has  the  following  analytical  solution: 

OX0,t)  I  rx 

\e/X\t>X() 

which  means  that  original  Eq.(4)  has  the  exact  analytical  solution 


*(*o>0  =  i 


n  +  2co 
en  ,t  < - -x, 


g  cos  6 
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2(0 


^ gcosO  0  ^  ^  2 CO 


g  cos  6 


-Xn 


i.e.  any  initial  condition  in  the  spatial  fetch  point  x0,  frequency  co  and  angle  0  is  growing 


exponentially  until  it  reaches  critical  value  e 


2  CO 

' - x 

geos# 


and  stays  equal  to  that  value  thereafter. 


(4) 


(5) 


(6) 


We  perfonned  numerical  tests  of  the  Eq.(4)  with  the  help  of  "rectangle"  numerical  scheme,  having 
second  order  approximation  in  space  and  time.  Such  scheme  is  unconditionally  stable,  and  the  only 
condition  of  its  correspondence  to  the  physical  reality  is  temporal  resolution  of  the  exponential  growth 

with  increment  y ,  i.e  time  step  of  the  numerical  integration  r  «  — .  This  condition  has  to  be  taken 

Y 

into  account  for  any  full-blown  numerical  simulation  of  Hassehnann  Eq.(l)  to  get  physically  sensible 
results. 


Fig. 3  shows  the  results  of  numerical  simulation  of  the  Eq.(4):  the  initial  condition  is  growing 
exponentially  until  it  reaches  critical  value,  and  stays  constant  afterwards. 
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Fig.  3:  Logarithm  of  3  spectral  energy  modes,  located  in  different  parts  of  the  fetch.  In  accordance 
with  analytical  solution  (6),  they  grow  exponentially  until  reach  the  critical  value. 

Fig.4  shows  total  energy  evolution  as  a  function  of  the  fetch  and  time  and  has  important  interpretation: 
negligence  of  the  nonlinear  interaction  term  leads  to  asymptotical  stationary  distribution  of  the  wave 
amplitudes  along  the  fetch.  However,  this  distribution  is  exponential,  and  such  distribution  of  energy 
has  never  been  observed  in  the  reality  (the  observed  is  the  power  one).  Therefore,  nonlinear  term  Snl  is 
at  least  absolutely  necessary  component  of  proper  physical  picture  of  the  ocean  wave  surface.  In  fact, 
Snl  is  the  dominating  source  tenn  in  Hassehnann  equation  with  respect  to  ;/ ,  as  it  was  shown  in  [1]. 


Fig.4:  Total  energy  as  a  function  of fetch  and  time 
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3.  Detection  of  swell  feedback  by  sea  background  through  theory  and  experimental  data 

The  global  visual  wave  observations  are  reanalyzed  within  the  theoretical  concept  of  self-similar  wind- 
driven  seas.  The  core  of  the  analysis  is  one-parameter  dependencies  of  wave  height  on  wave  period. 
Theoretically,  wind-driven  seas  are  governed  by  power-like  laws  with  exponents  close  to  Toba’s  one 
3/2,  while  the  corresponding  swell  exponent  (-1/2)  has  an  opposite  signature.  This  simple  criterion  was 
used  and  appeared  to  be  adequate  to  the  problem  of  swell  and  wind-driven  waves  discrimination.  This 
theoretically  based  discrimination  does  not  follow  exactly  the  Voluntary  Observing  Ship  (VOS)  data. 
This  important  issue  is  considered  both  in  the  context  of  methodology  of  obtaining  VOS  data  and 
within  the  physics  of  wind  waves.  The  results  are  detailed  for  global  estimates  and  for  analysis  of 
particular  areas  of  the  Pacific  Ocean.  Prospects  of  further  studies  are  discussed.  In  particular,  satellite 
data  are  seen  to  be  promising  for  tracking  ocean  swell  and  for  studies  of  physical  mechanisms  of  its 
evolution. 

Fig. 5  gives  a  graphical  summary  of  four  reference  cases  of  self-similar  evolution  of  wind-driven 
waves.  These  cases  are  shown  as  different  R,  tangents  of  one-parametric  dependencies 

H~7* 

height-to-period  in  logarithmic  axes.  Reference  cases  of  growing  wind  sea  are  shown  as  the  young  sea 
growth  at  pennanent  wave  momentum  production  (exponent  R  =  5/3  by  Hassehnann  et  al.  1976  [5]), 
growing  Toba’s  sea  (R  =  3/2 )  and  old  premature  sea  by  Zakharov  and  Zaslavsky  1983  [6]  with  R  = 

4/3. 


in  Ts 


Fig.  5  Reference  cases  of  wave  growth  as  one-parametric  dependencies  HS(TS).  Cases  of  Toba 
[1972]  R  =  3/2  law  and  swell  with  R  =  -1/2  are  shown  by  bold  lines.  Domain  where  the  asymptotic 
scheme  is  formally  invalid  (Badulin  et  al.,  2007  [7])  is  shaded. 
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Fig.  6  Dependencies  H(T)  and  their  power  law  fits  for  the  whole  World  Ocean,  1970-2007 for  (top) 
wind  waves  and  (bottom)  swell.  Lines  marked  as  A,  B,  C,  D  show  reference  power  laws  of  Fig.  5.  The 
exponents  of  the  experimental  fits  R  =  0.96  for  wind  sea  and  R  =  0.54  for  swell  are  found  to  be  quite 
far  from  the  reference  cases.  Totally,  36,356,695  reports  have  been  used  for  wind  waves  and 

31,041,169 for  swell  observations. 

The  height-to-period  dependencies  H(T)  derived  from  the  VOS  data  are  shown  in  Fig.6  for  visually 
delineated  wind  waves  (upper  panel)  and  swell  (bottom  panel).  The  difference  of  these  two  cases  is 
clearly  seen  in  tenns  of  exponents  R.  Thorough  analysis  of  the  experimental  data  uncovers  more 
physically  significant  difference  of  two  sea  waves  extremes.  Selecting  waves  in  wave  ages  and,  what  is 
more  important,  in  wave  periods  we  found  definite  indications  on  pumping  of  swell.  Results  of  such 
selection  is  illustrated  by  Fig.  7  where  histograms  of  exponents  R  estimated  for  20°  by  20°  coordinate 
boxes  are  given.  While  for  the  wind  waves  the  histogram  is  localized  near  a  value  more  than  1 ,  for 
swell  the  corresponding  distribution  is  quite  large  that  implies  a  variety  of  physical  mechanisms 
responsible  for  the  swell  dynamics. 
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Fig.  7  Histograms  of  exponents  R  of power  law  fits  of  H(T)  dependencies  calculated  monthly  for 
2 Ox  20  boxes  of  the  World  Ocean  (1481  dependencies  of  total  1 782  =  8  xl8xl2  coordinate  boxes)  for 
special  ranges  of  wave  periods:  (bottom)  wind  sea  with  T  =  5-10  s  and  (top)  swell  with 

periods  T  =  10-20  s. 

With  the  exponent  R  as  indicator  of  sea  wave  dynamics  we  make  a  conceptual  step:  we  study  a  link  of 
wave  heights  H  and  periods  T  rather  than  features  of  the  independent  data  sets.  The  separate  analysis 
basing  on  VOS  (Gulev  et  ah,  2004  [8])  or  satellite  data  (e.g.,  Zieger,  2010  [9])  gives  valuable 
information  on  ranges  of  wave  parameters  and  their  geographical  variability,  but  propose  quite 
primitive  vision  of  wave  dynamics.  Recent  attempts  to  combine  satellite  altimeter  observations  of 
wave  heights  and  mathematical  modeling  of  wave  dynamics  (Laugel  et  al.,  2012  [10]  )  propose 
reconstructions  of  full  spatiotemporal  structure  of  wind  wavefield.  This  study  is  based  on  extensive 
simulations  and  requires  thorough  theoretical  analysis.  The  interpretation  of  its  results  in  the  context  of 
burning  problems  of  sea  wave  physics  shows  a  good  prospect  for  further  study. 

4.  Analytical  and  numerical  proof  of  non-integrability  of  2D  free-surface  hydrodynamics 

We  studied  the  problem  of  integrabily  of  2D  hydrodynamics  of  fluid  with  free  surface  in  a  gravity 
field.  This  conjecture  was  formulated  in  1994.  Here  we  studied  the  integrability  for  potential  motion  in 
the  framework  of  Hamiltonian  truncated  equation  up  to  the  fourth  order: 
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H  =  - 
2 


\\SV2  +  'J/k'Jfdx  -  ^{(ky/f  -  (y/x  f)ijdx  + {y/  jfky/  +  y/k(t]k(>)ky/))}dx  + ...  (7) 


here  r/(x,t )  is  the  shape  of  a  surface,  y/(x,t )  is  a  potential  function  of  the  flow  given  at  the  surface  and 
g  is  gravitational  acceleration.  To  simplify  the  Hamiltonian  we  applied  canonical  transformation 
which  excludes  all  nonresonant  temrs  and  instead  of  (7)  we  deal  with  the  equivalent  Hamiltonian 


[h‘o),bdx  +  —  f 

■8b: 

T  •  f 

i 

J  k  2 J 

dx 

A 

db  *  ,  .  8b 

< - b  — 

ox  ox 


—  K\  b\ 


dx. 


(8) 


here  k  is  the  operator  multiplying  Fourier  harmonics  by  mudulus  of  wave  nambers.  Corresponding 
equation  of  motion  is: 


.8b  ,  /  p, 

l  —  =  (d,b  +  —  P 
dt  k  4 


b*-(b,2)-—(b*’—b2) 

/-s  v  /  ^  ’ 

ox  ox  ox 


-If* 

2 


b-K(\  b’  \2)--^{b'K(\b\2)) 
ox 


(9) 


This  equation  has  localized  breather-type  solution 

b(  x,  t)  =  B(x-  Vt)el(k()X~%t] ,  (10) 

where  k0  is  the  wavenumber  of  the  carrier  wave,  V  is  the  group  velocity  and  a>0  is  the  frequency 
close  to  ook  .  In  the  Fourier  space  breather  can  be  written  as  follow: 


Kit) 


-i{£lt+Vk) 


’4, 


(ii) 


00,. 


where  O  is  close  to  ^  •  This  solution  is  stable  and  does  not  radiate.  In  the  integrable  systems 

collisions  of  such  beaters  must  be  ellastic.  We  perfonned  numerical  simulation  of  collisions  of  two 
breathers  and  have  found  that  it  is  not  pure  ellastic,  see  Fig. 8. 


time=0.000000 


X 


Fig.  8  Surface  profile  of  two  breathers. 
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Also  we  studied  analytically  coefficient  of  6-waves  interaction  as  a  superposition  of  4-waves 
interaction.  For  integrable  system  it  must  be  equal  to  zero  on  the  resonance  manifold.  However  it  was 
found  that  it  does  not  vanish.  So,  both  numerical  and  anaylical  study  allow  us  to  conclude  that  2-D  free 
surface  hydrodynamics  is  not  integrable  system. 

5.  Dissipation  of  energy  during  the  wave-breaking  event 

One  of  the  most  challenging  tasks  in  improving  the  wave  forecasting  models  is  to  develop  an 
appropriate  dissipation  source  term.  The  wave  breaking  phenomenon  is  the  most  important  source  of 
dissipation  of  energy  for  waves.  Recently  if  has  been  found  that  the  modulational  instability  is  in 
principle  a  good  candidate  to  explain  the  wave  breaking  mechanism  in  deep  and  intennediate  water 
depth.  The  main  idea  is  that  a  wave  packet,  if  sufficiently  steep,  may  go  through  a  modulational 
instability  process  by  which  a  single  wave  grows  at  the  expense  of  the  other  waves.  If  the  limiting 
steepness  is  reached,  then  a  wave  breaking  is  observed.  Normally,  the  modulational  instability  is 
studied,  at  the  leading  order,  by  the  Nonlinear  Schrodinger  equation  and  its  exact  solutions.  Numerical 
simulations  of  the  fully  nonlinear  Euler  equation  have  also  been  perfonned  in  the  past;  however,  the 
major  limitation  of  these  approaches  is  that  they  fail  to  capture  the  final  stages  of  the  wave  breaking 
mechanism,  therefore  they  are  not  capable  of  furnishing  any  result  which  concerns  the  dissipation  of 
energy. 

In  the  present  work  we  have  performed  for  the  first  time  a  direct  numerical  computation  of  the  full 
Navier-Stokes  equation  for  two-phase  flow  (air  above  and  water  below  the  surface)  in  order  to  study 
the  wave  breaking  that  results  from  the  modulational  instability  mechanism.  If  the  steepness  of  the 
initial  wave  is  large  enough,  we  observe  a  wave  breaking,  as  shown  in  Fig. 9,  and  the  formation  of  large 
scale  dipole  structures  in  the  air,  see  Fig.  10.  Because  of  the  multiple  steepening  and  breaking  of  the 
waves  under  unstable  wave  packets,  a  train  of  dipoles  is  released  and  propagate  in  the  atmosphere  at  a 
height  comparable  with  the  wave  length.  The  amount  of  energy  dissipated  by  the  breaker  in  water  and 
air  is  considered  and,  contrary  to  expectations,  we  observe  that  the  energy  dissipation  in  the  air  is 
larger  than  the  one  in  the  water,  see  Fig.  1 1 . 

^  0.05 
6 

0.1 

Fig.  9  wave  breaking  during  the  modulational  instability  process 

0.1 

|  0.05 


Fig.  10  Dipole  formation  during  wave  breaking 
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Fig.  11  Dissipation  of  energy  during  multiple  breaking  events.  Dashed  line  corresponds  to 
dissipation  in  air  and  solid  line  to  dissipation  in  water. 

This  result  is  somehow  counter  intuitive:  most  of  the  experimental  research  so  far  have  measured 
dissipation  due  to  wave  breaking  by  looking  at  the  dissipation  of  kinetic  energy  in  the  water.  Possible 
consequences  on  the  wave  modeling  of  dissipation  source  tenn  and  on  the  exchange  of  aerosols  and 
gases  with  the  atmosphere  are  under  investigation. 

IMPACT/APPLICATIONS 

•  Confirmation  of  the  correctness  of  proposed  new  wind  input  term  for  limited  case  situation 
should  improve  operational  forecasting  models 

•  New  analytic  solution  for  Hassehnann  equation  without  nonlinear  tenn  will  help  to  test 
operational  models  against  correctness  of  the  physics  of  described  phenomena 

•  Detection  of  swell  feedback  by  sea  background  through  theory  and  experimental  data  helps  to 
understand  complex  nonlinear  mechanisms  of  different  scales  surface  waves  interaction  in  the 
ocean 

•  Analytical  and  numerical  proof  of  non-integrability  of  2D  free-surface  hydrodynamics  will  help 
to  build  simplified  and  fast  nonlinear  models  of  surface  waves 

•  Discovery  of  large-scale  dipole  structures  in  the  air  appearing  due  to  wave-breaking  might  alter 
current  presentation  of  dissipation  source  terms  in  operational  models 

RELATED  PROJECTS 

None 
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